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We derive a fully relativistic, self-similar solution to describe the propagation of a 
shock along an exponentially decreasing atmosphere, in the limit of very large Lorentz 
factor. We solve the problem in planar symmetry and compute the acceleration of 
the shock in terms of the density gradient crossed during its evolution. We apply our 
fT} ■ solution to the acceleration of shocks within the atmosphere of a HyperNova, and show 

that velocities consistent with the requirements of GRB models can be achieved with 
exponential atmospheres spanning a wide density range. 
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*-«^ . 1. Introduction 

In a SuperNova (SN) explosion, the shock wave must ultimately emerge from the body of the 
star, and begin to propagate down the exponential density gradient of the stellar atmosphere. In 
Newtonian fluid dynamics, the propagation of such a shock is described by a self-similar solution 
(Raizer 1964; Grover and Hardy 1966; Hayes 1968). This self-similar solution is unavoidable 
(Raizer 1964). In fact, it turns out that, as the shock accelerates, a sonic point is formed, 
separating matter located immediately behind the shock from the flow's initial conditions; this 
prevents pressure waves (i.e., causal information) to reach post-shock material from the the area 
where the flow's initial conditions are set. Thus all flows will converge to the same solution, 
irrespective of their initial conditions. The lack of dependence upon initial conditions severely 
constrains the ensuing flow, by restricting the number of parameters upon which the shock 
evolution may depend; in fact, this constraint is so strong that only a single (self-similar) solution 
exists. 

In the Newtonian solution, it is found that the shock velocity increases exponentially, in such 
a way that the shock reaches spatial infinity in a finite time. This is clearly impossible when 
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account is taken of Special Relativity, and therefore the Newtonian solution must break down at 
high shock velocities, as expected. The question of precisely how the shock evolves in relativistic 
conditions has not been investigated so far. 

This problem has acquired a new urgency within the HyperNova model (Paczyhski 1998; 
MacFadyen & Woosley 1999) for Gamma Ray Bursts (GRBs). In fact, Meszaros and Rees (2001) 
have shown that, despite the large assumed energy release of the central engine, the outwardly 
moving shock in the large star hypothesized to give rise to the GRB can only reach a Lorentz factor 
of Fi « 10 at the end of the H-envelope, and they had to invoke (without explicit computations) 
shock acceleration down the star exponential atmosphere to reach the required Lorentz factors of 
Tf m 100 — 300. The existence of a self-similar temporal structure in GRBs is also suggested by 
an analysis of their power spectra (Beloborodov, Stern, & Svensson 2000; see also Sivron 1998). 

It is the purpose of this paper to derive a fully (Special) relativistic self-similar solution for 
the problem of shock propagation in an exponential atmosphere (Section 2), and to discuss its use 
in the HyperNova model for GRBs (Section 3). 



2. Relativistic self similar flow 

We shall assume that a cold (T = 0) material is stratified with density distribution 
p = p ex.p(—kx); the shock is supposed to move toward x = +oo, thus with positive velocity 
v > 0. The symmetry is assumed to be planar, since both the length scale of stellar atmospheres 
and the total extension of the atmosphere are much smaller than the stellar radius. In the 
Newtonian analog, the shock speed is set by a purely dimensional argument (e.g. Chevalier 1990): 
V = a/kt, with a an adimensional constant to be determined. In the special relativistic problem, 
dimensional arguments alone fail because the presence of the light speed c allows the construction 
of a new adimensional quantity (kct), thusly spoiling the above argument. We can however recover 
from this impasse by appealing to both dimensional and covariance arguments. The time-like 
part of the shock four-speed is of course determined by the identity UpJJ^ = —1, while the 
space-like part must be built in covariant fashion from the available quantities. Defining a four 
vector k^ = (0, k, 0, 0), we can then write for the spatial part of U 

dX a ak a 

U a = — = (1) 

ds k^s v ; 

where k^ = k 2 > 0. This is the only sensible solution to our dimensional/covariant problem. We 
see from this that, as long as we use proper time, the structure of the problem is identical to that 
of the Newtonian analog. In particular, the shock reaches spatial infinity within a finite proper 
time. We may choose s = for the moment when this occurs, so that the flow is restricted to 
s < 0. Incidentally, note that this implies a < 0. Physically, this makes perfect sense: as the shock 
accelerates, its proper time is contracted by its Lorentz factor (r) with respect to the fluid time. 
The shock speed in terms of the fluid time can be obtained by remembering that dX/ds = vT, and 
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substituing in the above equation one finds (c = 1 from now on) 



Taking the time derivative of both sides we obtain 



a \dv vT 
which can be immediately integrated to yield 

kt , fl + v\ l l 2 1 

= log h constant . (4) 

—a \1 — v J v 

The low speed limit v <C 1 is kt/a « in agreement with the Newtonian solution. The 
hyperrelativistic limit (v — ► 1) is 

^-^logr/r, (5) 

—a 

where we have introduced an initial shock Lorentz factor Tj at time ti, for ease of use in the 
future. Since the shock moves essentially at speed 1, we can rewrite the lhs of the above equation 
as log(p/p i ) 1 / a , from which we see that the initial and the final (i.e., when the shock leaves the 
exponential atmosphere) shock Lorentz factors are related by 

The above clearly shows that all we have left to do is to determine the value of the parameter 
a. To this purpose, we must consider the fluid equations for the post-shock material, which we 
have found it convenient to take in the form given by Blandford and McKee (1976, their Eqs. 
14-15) 

^0g(eV) = -4| (8) 
dn d . 

a+ s ("0 = o (9) 

where d/dt = d/dt + vd/dx is the convective derivative, and we have dropped curvature terms 
(which are contained in the Blandford and McKee equations), to keep with our planar symmetry 
approach. Here e is the local energy density, 7 is the local fluid Lorentz factor, as seen from the 
reference frame of the unshocked material (to be distinguished from T, the shock Lorentz factor 
in the same reference frame), and n is the baryon number density always in the reference frame 
of the unshocked fluid. The above equations assume an hyperrelativistic equation of state for the 
post-shock material of the form p = e/3, which is correct in the limit V — ► 00. Indeed, for the 
large shock Lorentz factors appropriate to this problem (r ^ 10) Heavens & Drury (1988; see 
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also Iwamoto 1989) showed that the adiabatic index of the flow is within 1% of its asymptotic 
relativistic value of 4/3, irrespective of the upstream temperature and plasma composition (see 
especially their Fig. 2, and Eqs. 25-26). It should be remarked here that Eq. (JgJ) remains valid 
whatever we assume for the equation of state, but the precise value of a depends instead on the 
particular choice of the equation of state. 

The boundary conditions for this problem are provided by Taub's jump conditions (Taub 
1948), which again we take in the hyperrelativistic limit given by Blandford and McKee (1976): 

e 2 = 2rV , 7 | = ^r 2 , n 2 = 2rV . (10) 

Here the subscripts 1,2 refer to pre- and post- shock quantities, respectively; n and 7 are 
always defined in the unshocked fluid frame, and e in the comoving frame. The enthalpy of the 
unshocked material, wi, equals e\ since this material is assumed to be cold. In our problem, e\ is 
not a constant, because the atmosphere is stratified. We have 

wi = ei = pi = P o exp{-kX) = Po{T/Ti) a (11) 

so that, ultimately, 

62 = 2 (r?) r2+ ° = 2q ° r2+a (12) 

which is the form we shall use in the following. An identical argument shows that 

n 2 = 2^T 2+a = 2z T 2+a , (13) 

i 

where n D = p /m. 

In order to search for a self-similar solution, we need to assume a form for the similarity 
variable. In the Newtonian analog, this is clearly £ = k(x — X), where X is the instantaneous 
shock location. In this problem, we take 

£ = k(x- X)T 2 . (14) 

The rationale for this is that, since post-shock material has a Lorentz factor which is v2 times 
smaller than the shock's, the post-shock material falls behind the shock by an amount oc 1/T 2 
as the shock covers an exponential length. Please notice that, with our notation, £ < for the 
shocked material. The post-shock quantities 7 2 , e and n can then be taken of the form 
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g(0T 2 , e = q o R(0r 2+a , n = z o T 2+a N(0 (15) 



with boundary conditions given by Eqs. ([!(]) (the second one), (|i~2"D, (O) as 



g(0) = ~ , R(0) = 2 , N(0) = 2 . (16) 
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Substitution of Eqs. (0), (||) into Eqs. (@), (|) and (|) yields, aft er some algebra (the dot 
indicates derivation with respect to £) 

A = -g(2/a + 1)(1 - 4g/q + 1/g) + (10/q + 3) 

12 3/2-3/2 5 -6£/a + (l-g/2 + 2£ 5 /a)(l-4£/a + l/g) 1 ' 

4 = |(l/5- 1/2 + 2£/a) + (2/a + l) (18) 

5 12 

_ = i (IQ) 

N g(l - 4£/a) - 1 ' V j 

From these, we see what fixes a: the denominator of the rhs of Eq. (|l~7| ) goes to zero at a critical 
point £ c and, unless the numerator simultaneously does the same (which will only occur for a 
special value of a), a non-integrable singularity will ensue. The same phenomenon occurs in the 
Newtonian analog, where it has been shown that this critical point is a sonic point. It is exactly 
because there is a sonic point that a self-similar solution can develop: in fact, material between 
the shock and the sonic point is not in causal contact with post-sonic point material, and thus its 
properties cannot be determined by the problem's initial conditions. 

To keep the problem non-singular, we impose that the numerator and denominator of Eq. 
( JT7D vanish simultaneously, obtaining 

g ^ c )(l-^ c /a) = ^±^ (20) 
2 + a 

from the numerator, and 

5 fe)(l-4£ c /a) = i^? (21) 

for the denominator. The conditions above lead to two solutions for a, one positive and another 
negative. The positive solution must be discarded because otherwise the shock would not reach 
spatial infinity for s — > (Eq. [l]). We are therefore left with only one sensible solution, 

— 12 — \/l92 

a = « -4.309401 . (22) 

6 

The location of the sonic point, £ c , where both the numerator and the denominator vanish, is 



determined through numerical integration of Eqs. (|17|), ( |i8[ ) and (19). This yields £ c w —0.46. 
Also note that g(£)(l — 4£/a) < (?(0) = 0.5, and therefore the denominator of Eq. ( |i~9| ) is always 
well-behaved. The full numerical solution for <?(£), R(0 and iV(£) is illustrated in Fig. 1. 



3. Application to GRBs 



It can be seen from Eq. (0) that the necessary asymptotic Lorentz factor required for a proper 
modelling of the properties of GRBs, i.e. ^ 150 (remember that the Lorentz factor of the 
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matter is only l/y/2 that of the shock) can be reached by crossing a modest factor of 10 5 density 
decrease in the exponential atmosphere (assuming Tj ~ 10). Yet, the total energy of material 
moving at these large Lorentz factors is only modest. In fact, let us compute the distribution of 
kinetic energy with baryon number. Consider a cylindrical fluid element with surface area A (with 
the normal along the direction of fluid motion x), and length dx along the direction x. This element 
will have a bulk kinetic energy of dE = Amc 2 n(x)-y(x)dx = Ap c 2 T 1+a T^ a N(^)g 1 ^' 2 (^)k~ 1 d^. 
What we want to know is the energy distribution at the moment in which the shock gets out of the 
exponential atmosphere. Let us indicate with Tf the shock Lorentz factor at that moment. Then 

* 2 ^-'r-«rt+«^^ (23) 

dy 1 1 g 

which can be coupled to the solution g = to yield a parametric representation of jdE/dy, the 
distribution of kinetic energy with respect to the final Lorentz factor. The adimensional part of 
the function in the previous equation is plotted in Fig. 2. The numerical factor can be estimated 
using Eq. (0): 

7 ^ = io« erg M^-V ( J ^~) 1 (~Y (^-Y +1 2iV(0g3/2(g) (24) 

1 dy Hl0 13 cnj ^10 12 cmy U0~ 9 cm" 3 / \rj \15qJ g ' 1 ' 

In this expression, the values for ru (the edge of the H-shell), pi (the matter density at the end 
of the H-shell, and thus presumably at the beginning of the exponential atmosphere), and r, (the 
shock Lorentz factor at the end of the H-shell = beginning of the exponential atmosphere) are 
taken from Meszaros and Rees (2001). Taking the adimensional factor from Fig. 2, we see that 
the kinetic energy falls short of the average isotropic burst energy (Schmidt 1999) by about three 
orders of magnitude (note that here what we are computing is the isotropic energy). 

Before showing how to come out of this impasse, we need to consider which part of Fig. 2 
can actually be obtained in a realistic model. The reason for this limitation comes from the fact 
that the idealized model presented here has been propagating down an exponential atmosphere 
forever, while in a realistic model only a finite amount of matter can have converged onto this 
self-similar solution, given the finite dimension of the star. For this reason, we have plotted in Fig. 
2 three ticks, which correspond to the present position of matter which was located, before being 
reached by the shock, at 5, 10, 15 exponential scale-lengths from the shock's present location. If 
we think that the initial atmosphere extends for 5, 10, 15 exponential scale-lengths, then we can 
believe the part of Fig. 2 located to the right of their respective tickmarks. Leftward of them, the 
true, physical solution will depart from the one shown here, and the total kinetic energy at the 
corresponding values of Lorentz factor will be much smaller than the values plotted (obviously, 
since the physical solution has finite mass and energy, which is not true for the idealized one). 

Let us now suppose that the putative burst progenitor possesses a wide exponential 
atmosphere, spanning some ^ 9 orders of magnitude in density. Then, for an initial Lorentz factor 
Ti ~ 10, the shock Lorentz factor at the end of the atmosphere is ~ 1000 (Eq. |6|); once again, the 
overall energy of matter moving at this speed would be modest (~ 10 46 erg), but Fig. 2 shows that 
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most of the energy (before the cutoff implied by the finite extent of the atmosphere) will come 
out at a Lorentz factor ~ 0.2 of the shock's factor, i.e. at ^matter ~ 200. With the values above, 
the total energy then amounts to ~ 6 x 10 52 erg, in reasonable agreement with Schmidt's (1999) 
isotropic estimates^. It should be noticed that inspection of the numerical solution for large values 
of the distance parameter £ shows the material to be cold (e oc n), so that there will be negligible 
further acceleration of these slower shells by pdV work, and r mat t er ~ 200 remains a good estimate 
of the coasting Lorentz factor. 

To summarize, the application of our solution to shock acceleration in the atmosphere of 
a massive star has shown that, in atmospheres with small density range, the amount of energy 
carried by material accelerated to the typical GRB Lorentz factors falls short of the GRB required 
energetics. However, for stars with a wide density range in their atmospheres (;> 10 9 orders 
of magnitude), a sufficient quantity of energy is carried by later shells of material moving at 
the typical GRB Lorentz factors. In this model, the early emission would be dominated by an 
ultra-hard component, due to the very large-T shells slowing down in the ISM. This is at least 
qualitatively consistent with virtually all bursts studied in some detail by BeppoSAX (Frontera et 
al., 2000, see especially their Fig. 2), which exhibit, in their first few seconds only, spectra peaking 
beyond the instrument's (GRBM) observing limit of 700 keV. 

RP thanks the Osservatorio Astronomico di Roma for its kind hospitality during the time 
that this work was carried out. 
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Fig. 1. — Numerical solution of Eq. (|l7|) , (|l8|), and (|i~9l) with a given by Eq. (p^). 
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Fig. 2. — The quantity 2A r (£)g(£) 3/ ' 2 /g (the adimensional part of jdE/dy, Eq. |23| ), plotted as a 
function of 1//2 . The ticks in boldface indicate the present position of matter which, before 
being shocked, was located at 5, 10 and 15 exponential scale lengths k^ 1 . 



